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Abstract 



We consider the problem of pricing a contingent claim whose payoff depends on several 
sources of uncertainty. Using classical assumptions from the Arbitrage Pricing Theory, the 
theoretical price can be computed as the discounted expected value of future cash flows 
under the modified risk-neutral information process. Although analytical solutions have 
been developed in the literature for a few particular option pricing problems, computing the 
arbitrage prices of securities under several sources of uncertainty is still an open problem in 
many instances. In this paper, we present efficient numerical techniques based upon Monte 
Carlo simulation for pricing European contingent claims depending on an arbitrary number 
of risk sources. We introduce in particular the method of quadratic resampling (QR), a new 
powerful error reduction technique for Monte Carlo simulation. Quadratic resampling can be 
efficiently combined with classical variance reduction methods such as importance sampling 
to further improve the accuracy of the estimate. Our numerical experiments show that the 
method is practical for pricing claims depending on up to one hundred underlying assets. 
We also describe an implementation of the method on a massively parallel supercomputer, 
yielding two orders of magnitude of performance improvement over the same implementation 
on a desktop workstation. 



Resume 



Nous etudions le probleme de 1'evaluation d'un actif conditionnel dont la valeur depend de 
plusieurs sources de risque. Avec les hypotheses classiques de la theorie de 1' arbitrage, le prix 
theorique peut etre calcule comme l'esperance mathematique actualisee des flux financiers 
futurs pour le processus modifie neutre vis-a-vis du risque. Bien que des solutions analytiques 
existent pour certains problemes particuliers d'evaluation d'options, 1'evaluation numerique 
du prix d' arbitrage des actif s multidimensionnels reste un probleme ouvert a ce jour dans de 
nombreux cas. Nous proposons des techniques numeriques efficaces basees sur la methode de 
Monte Carlo pour calculer le prix d' actif s Europeens dependants d'un nombre arbitraire de 
sources de risque. Nous introduisons en particulier la methode du reechantillonnage quadratique 
(RQ), nouvelle technique puissante de reduction d' erreur pour la simulation de Monte Carlo. Le 
reechantillonnage quadratique peut etre combine avec les methodes classiques de reduction de 
variance telles que l'echantillonnage par importance pour reduire encore l'erreur d'estimation. 
Les resultats experimentaux montrent l'efficacite de cette approche pour evaluer des actifs 
conditionnels dependants d'une centaine d'actifs sous-jacents. Nous presentons aussi une 
implantation de la methode sur une machine massivement parallele, permettant de gagner deux 
ordres de grandeur en vitesse de calcul par rapport a une implantation sur station de travail. 
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1 Introduction 

1.1 Background 

Since the seminal work of Black and Scholes [Black and Scholes, 73] and Merton [Merton, 73] 
in the early 1970s, the arbitrage principle underlying option valuation theory has been extended 
to a broad range of other financial instruments (see e.g. [Ross, 76], [Cox and Rubinstein, 85]). 
Indeed, any security whose returns are contractually related to the returns on some other 
security or group of securities can theoretically be valuated using the same arbitrage principle. 
This is the case in particular of warrants, convertible bonds, but also common stocks, 
ordinary bonds, and most other types of contractual instruments. In some cases, explicit 
closed form analytical formulas for the computation of the arbitrage price can be derived 
from this theory. In particular, the original paper of Black and Scholes provides a closed 
form solution for a European option on a single common stock. Unfortunately, few other 
cases can be solved analytically, and computing the arbitrage price often requires numerical 
simulations. Following an idea initially presented in an early edition of [Sharpe, 85], Cox, 
Ross, and Rubinstein [Cox, Ross, and Rubinstein, 79] [Cox and Rubinstein, 85] developed a 
discrete model for the valuation of an American option on a single stock that can be easily 
computed numerically. However, the effective implementation of the arbitrage principle is 
not always such an easy task, and may sometimes become intractable. For example, when a 
contingent claim depends on many sources of uncertainty, there is no known general tractable 
algorithm for computing accurately the arbitrage price. 

1.2 Multidimensional pricing models 

There are several reasons motivating the development of efficient methods for multidimensional 
contingent claim pricing. We briefly review some of them below, but the list is certainly not 
exhaustive. The tremendous development of financial engineering during the past decade can 
be expected to continue, and new types of securities requiring multidimensional modeling are 
likely to appear at a sustained pace in the future. 

• Complex Over The Counter warrants 

Many instruments depending on several underlying securities are being proposed on 
Over The Counter markets. One can buy options on the best performance of several 
securities, or on the spread between different market indexes. In foreign exchange 
markets, complex instruments are being proposed that enable investors to diversify their 
portfolios internationally while protecting them from the relative variations of several 
different currencies. In commodity markets, long-dated options on portfolios of futures 
contracts with various maturities are being proposed to protect investors from long term 
evolutions of the underlying commodity and of its convenience yield. 

• Path-dependent instruments 
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Specialized institutions are now proposing options that depend not only on the current 
value of the underlying asset, but on the whole history of underlying asset prices 
until expiration date. Well known examples are lookback options (e.g. options on the 
minimum value of the underlying for a given period of time), and Asian options (e.g. 
options on the average value of the underlying for a given period of time). These 
path-dependent instruments can be viewed as multidimensional contingent claims on 
the sequence of underlying asset prices until the expiration date. Other examples of 
path-dependent instruments are Mortgage-Backed Securities (MBS), which depend on 
the past history of interest rates, and certain life insurance contracts, such as Single 
Premium Deferred Annuities (SPDA) and Guaranteed Interest Contracts (GIC). 

• Interest rate related instruments 

A new generation of term structure models [Heath, Jarrow, and Morton, 92] has recently 
been developed. These models can account for non-parallel evolutions of the term 
structure, i.e. discrepancies in the evolutions of short and long-term rates. According to 
these models, the whole term structure of interest rate is a multidimensional underlying 
asset for interest rate contingent claims. Examples of interest rate contingent claims 
include all kinds of fixed-income securities such as bonds and convertible bonds, bond 
options, swaps, caps, floors, collars, swaptions, but also mortgage-backed securities and 
SPDAs [Fabozzi, 87]. 

• Models with stochastic parameters 

Is is well known that the volatility of the underlying asset cannot be assumed constant 
for long-term options. Recently, new option pricing models have been developed that 
attempt to take into account the stochastic nature of volatility (see e.g. [Wiggins, 87], 
[Dothan, 87], [Hull and White, 88]). In these models, the option price depends on two 
sources of uncertainty: the underlying asset itself, and its volatility. Similarly, the 
interest rates cannot be assumed constant for a long period of time. Therefore, pricing 
long-dated options requires modeling the uncertainty arising for interest rate variations. 
Examples of commodity or stock option pricing models accounting for uncertainties 
in stochastic interest rates can be found in [Merton, 73], [Jarrow, 87]. Of course, as 
described above, fixed income securities and related derivative instruments also require 
modeling of the stochastic nature of interest rates. More generally, any parameter in the 
stochastic model of the underlying asset that cannot be assumed predictable constitutes 
a new source of uncertainty that cannot be neglected for long-dated instruments. 

• Quality delivery options in futures contracts 

The seller of a futures contract must deliver at expiration date some predefined underlying 
asset. Sometimes, the contract specifies a series of acceptable underlying assets, in 
which case the seller will deliver the cheapest one. In essence, the futures contract is 
a multidimensional option on the minimum value of all the possible deliverable assets. 
Consideration for the quality option is especially important in the case of interest rate 
futures, where a number of bonds with different durations and characteristics can often 
be delivered [Cheng, 87], [Boyle, 89]. 
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• Pricing of insurance contracts 

The Arbitrage Pricing Theory can be applied to the pricing of insurance policies. In 
particular, as described above, some life insurance policies can be viewed as pure interest 
rate contingent claims (see e.g. [Stanton, 89]). More generally, option pricing theory 
can be applied to the valuation of property /liability insurance contracts, where the strike 
price of the option corresponds to the deductible of the insurance policy. There is a 
growing academic literature on these subjects, that can be traced back to the work of 
[Merton, 77] and [Smith, 79]. Other applications are presented in [Kraus and Ross, 82], 
[Doherty and Garven, 86], [Cummins, 88] [Shimko, 92]. 

• Assets and Liabilities Management 

Large corporations bear many different kinds of risks in the course of their activities. 
Some of these risks can be hedged on financial markets. In particular, the currency 
exchange risks and the interest rate risks can be hedged in many cases. For financial 
institutions such as major banks and life insurance companies, this financial risk is 
the only major source of uncertainty of the business. In order to minimize this risk, 
the corporation must immunize itself as much as possible against variations of the 
economic value of the firm as a function of the sources of uncertainty. This economic 
value is simply the difference between the present value of assets and the present 
value of liabilities. Given a model of the commercial and financial policy of the firm, 
these present values can be computed as contingent claims on the economic variables 
representing the risk (e.g. interest rates and currency exchange rates). 

• Capital budgeting models in corporate investment decision making 

When assessing the opportunity to invest in a new project, a financial manager must 
determine the Net Present Value (NPV) of the project, i.e. the economic value of all the 
possible future cash flows generated by the project. Several economic and commercial 
variables are generally relevant for building a model of the probability distribution of 
the future cash flows. In particular, when building the business plan for introducing a 
new product to the market, one must estimate the possible evolutions of the parameters 
associated with the product, as functions of the predefined marketing policy. These 
parameters can be the market size, the market share of the company, the fixed and 
variable costs, the expected price. Also, global economic variables such as gross 
national product, rate of inflation, or interest rates may be relevant. According to 
[Brealey and Myers, 91], the net present value can be computed as the expected value of 
all future cash flows, discounted at the opportunity cost of capital. Computing the NPV 
can be viewed as pricing a multidimensional contingent claim on the relevant economic 
and commercial parameters. See [Mason and Merton, 85] for a review of applications 
of option pricing theory to corporate finance. 

1 .3 European versus American instruments 

For pricing purposes, financial assets can be divided into two majors classes. The first class is 
that of assets whose future cash-flows cannot be influenced by decisions from the holder taken 
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after the purchase date. Most assets traded on financial markets belong to this class. We will 
call European instruments all financial assets belonging to this first class. In particular, stocks, 
bonds, futures contracts, European options, swaps, caps, floors, mortgage-backed securities are 
European instruments. The second class is that of assets whose cash flows can be influenced 
a posteriori by the holder. American options belong to this class. As another example, the 
crediting policy of a Life Insurance company selling SPDAs greatly influences the present 
value of the liabilities of the company. This crediting policy can be adjusted by the company 
after signature of the contracts with the policyholders. Therefore, the liability associated with 
the sale of SPDA can be viewed as an American security. We will call American instruments 
all financial assets belonging to this second class. 

Following the general theory of arbitrage pricing, the theoretical price of a European contingent 
claim is the discounted expected value of its future cash flows under the so-called "risk- 
neutral" probability distribution of the underlying economic factors [Harrison and Kreps, 79] 
[Harrison and Pliska, 81], [Duffie, 88], [Karatzas and Shreve, 88]. Mathematically, computing 
the arbitrage price reduces to computing an integral (sum) over the space of the underlying 
economic factors. When the dimension of the space of the underlying economic factors is 
small, standard techniques for numerical integration can be used. In some cases, the integral 
can even be computed analytically (e.g. Black-Scholes formula). However, the computational 
complexity of evaluating the integral is clearly exponential in the dimension of the space. 

In this paper, we present efficient numerical techniques based upon Monte Carlo simulation 
that enable us to compute accurate approximations of realistic European pricing problems in a 
time quadratic in the number of underlying economic factors. Combining classical techniques 
of importance sampling with a new approximation method called quadratic resampling, we 
can compute the theoretical prices of the most complex instruments within tens of seconds 
on a PC or a workstation, and within a fraction of a second on currently available massively 
parallel supercomputers. The method has been tested for pricing problems with up to 100 
underlying factors. 

The price of an American claim is the maximum over all possible cash flow monitoring strate- 
gies of the associated present values of cash flows. For example, the price of an American 
option is the maximum over all possible early exercise strategies of the corresponding present 
values. Since the space of cash flow monitoring strategies is generally huge, direct maximiza- 
tion of the present value is rarely practical (see [Bossaerts, 89] for a discussion). However, 
when the underlying economy is modeled as a Markov process, one can use the Bellman 
principle of dynamic programming to compute the optimal monitoring strategy. American 
options are typically priced using a discrete approximation of the dynamic programming 
principle. This is the case in particular of the CRR model [Cox, Ross, and Rubinstein, 79] 
for American stock option pricing. This approach becomes however impractical when the 
underlying economic space has many dimensions, since the dynamic programming algorithm 
requires a memory space exponential in the number of dimensions. This fact is known as the 
"curse of dimensionality" problem [Bellman, 57] for dynamic programming. Unfortunately, 
the techniques developed in this paper for European claim pricing do not generalize simply to 
the case of American claims. 
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This paper is organized as follows. In Section 2, we relate our contribution to previous work 
in Monte Carlo integration and asset pricing. In Section 3, we recall the usual assumptions 
on the stochastic processes governing the evolution of securities prices, and the main result 
of the Arbitrage Pricing Theory. In Section 4, we briefly review Monte Carlo techniques 
for the numerical valuation of multidimensional integrals, and describe in particular modern 
importance sampling methods. In Section 5, we present the method of quadratic resampling. 
In Section 6, we present numerical results demonstrating the efficiency of the quadratic 
resampling method when combined with appropriate importance sampling schemes. We also 
describe a massively parallel implementation of our Monte Carlo procedure. Finally, in 
Section 7, we discuss the capabilities and limitations of the approach. 

2 Relation to other work 

Good reviews of the Monte Carlo method and different variance reduction techniques such 
as antithetic variables, covariates, stratified sampling, importance sampling can be found 
in many sources such as [Hammersley and Handscomb, 64], [Zaremba, 68], [Haber, 70], 
[Kalos and Whitlock, 86], [Decker, 91] and references thereof. [Stroud, 71] is a comprehensive 
source of results regarding multidimensional Gaussian quadrature rules. All our experiments on 
Gaussian rules are based on methods proposed in [Stroud, 71]. The additive importance sampler 
presented in this paper is described in [Dantzig and Glynn, 89] and [Dantzig and Inf anger, 91] 
in the context of portfolio optimization and other operations research applications. 

The application of the Monte Carlo method to option pricing was first presented in [Boyle, 77], 
in the context of claims contingent to a single underlying asset. It has then been used by 
several authors for the valuation of path dependent contingent claims. In particular, the 
method has been used for pricing mortgage-backed securities (see [Schwartz and Torous, 89], 
[Hutchinson and Zenios, 91]). Interestingly, the model described in [Schwartz and Torous, 89] 
could be solved faster by direct finite-difference approximation of the arbitrage partial 
differential equation. Indeed, although the problem is initially path-dependent as a function 
of interest rates, it can be made path-independent by including two additional state variables 
in the model. These new variables aggregate all the necessary information on past values of 
interest rates. The final model is four-dimensional, which is still tractable for direct numerical 
integration. This technique of including new state variables to remove path dependence is 
described in much detail in [Stanton, 89]. Many problems of path-dependence can be solved 
that way. The Monte Carlo method, although much slower than direct numerical integration 
for problems with few dimensions, is very flexible and simple to implement. This is why 
it has been used in many problems of path-dependence that could be solved faster by direct 
numerical integration. Another application of the Monte Carlo method to the valuation of 
Single Premium Deferred Annuities contracts is described in [Stanton, 89]. 

The valuation of options depending of several underlying assets has been studied by several 
authors. [Brennan and Schwartz, 79] addresses the problem of pricing options under two 
sources of risk by direct finite-difference approximation of the generalized Black-Scholes 
equation. In this example the two sources of risk are the short term and the long term 
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interest rates. The approach is clearly limited to a few assets, since the memory space 
requirements and the computation time are both exponential in the number of underlying 
assets. [Boyle, Evnine, and Gibbs 89] developed a multinomial lattice method for pricing 
multidimensional options, in the spirit of the approach outlined in [Cox and Rubinstein, 85]. 
According to the authors, the computation becomes very burdensome for more than two assets. 
In fact, the multinomial lattice approach can be viewed as a finite-difference approximation 
of the generalized Black-Scholes equation using an explicit Euler scheme and an appropriate 
change of variables [Brennan and Schwartz, 78]. 

[Stulz, 82] presents an analytical solution to the problem of pricing a European option on 
the maximum or minimum of two underlying assets. The analytical solution is generalized 
in [Johnson, 87] to the case of an arbitrary number of assets, taking as given the cumulative 
multivariate normal distribution function. [Boyle, 89] and [Boyle and Tse, 90] developed an 
approximate method for the same problem. Although the problem is solved analytically 
in [Johnson, 87], the approximate method does not require preliminary computation of the 
cumulative multivariate normal distribution function. To the best of our knowledge, the above 
problem of option pricing on the maximum or minimum of several assets is the only case 
reported in the literature where solutions to high-dimensional (say, more than 5) problems 
have been developed. 

Our contribution to the problem of multidimensional contingent claim pricing reported in this 
paper is twofold: 

-First, we show that a proper implementation of the Monte Carlo method makes the problem 
of European contingent claims pricing tractable, even when the number of underlying assets is 
very large (up to 100 in our experiments). 

-Second, we present a set of three techniques for improving the accuracy of the Monte Carlo 
estimation. 

• We introduce an original error reduction technique, quadratic resampling, that consider- 
ably improves the accuracy of Monte Carlo integration. Quadratic resampling can be 
used for any multidimensional integration problem, and is computationally practical for 
integrals with up to a few hundred dimensions. It can be combined with classical vari- 
ance reduction techniques such as importance sampling to further improve the accuracy 
of the computation. It is shown that, given a random vector having arbitrarily many 
coordinates and an arbitrary joint distribution, the expected value of any polynomial of 
degree two or less in the vector coordinates is computed exactly when using quadratic 
resampling. 

• We introduce a deterministic Quasi-Monte Carlo method which does not require the use 
of pseudo-random numbers. For integrals in spaces with few dimensions, this Quasi- 
Monte Carlo sampling method is considerably more accurate than pseudo-random 
sampling. This method combines the advantages of direct numerical integration in 
low-dimensional spaces with those of Monte Carlo integration in high-dimensional 
spaces. 
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• We show that with currently available massively parallel technology, multidimensional 
contingent claim prices can be computed two orders of magnitude faster than on a desktop 
workstation. Therefore, since the accuracy of the Monte Carlo estimate is proportional 
to the square root of the number of samples, the massively parallel implementation is 
one order of magnitude more accurate than the workstation implementation for a given 
computing time. 

3 Arbitrage pricing of European securities 

The arbitrage pricing theory is described in many textbooks. We refer the reader to [Duffie, 92] 
for a comprehensive presentation. 

3.1 Diffusion model of information process 

We model the economy as a finite-dimensional vector of real-valued state variables = 
(xi(t), . . . ,x„(t)), called factors, representing all the information available to investors at 
time t. By definition, X(t) is known at time t. However, X(t) does not usually give enough 
information to investors to allow perfect forecasting of the future evolutions of security 
markets. Therefore, X(r), r > tis a stochastic process. 

By construction, X[t) is a Markov process, i.e. future values of vector X are independent of 
past values of vector X, conditionally to the knowledge of X today. Indeed, if future values of 
X were to depend not only on X(t) but also on past values, this would mean that X(t) does not 
represent all relevant information at time t, contradicting its definition. 

Since X(t) represents all information available to agents at time t, in a frictionless market, 
prices of securities must be deterministic functions of time and X(t). It is said that securities 
are contingent claims on the state variable X(t). 

For the sake of simplicity, we will assume that the information process X[t) is a diffusion 
process. However, our results on Monte Carlo valuation described in the next sections apply 
to more general types of stochastic processes. If X(t) is a diffusion process, it is a solution of 
a stochastic differential equation of the type: 

dX = A(X( t),t)dt + B(x(t),t)dW (l) 

or equivalently: 

n 

Vi G [l,n], dxi = a,(xi, . ..,*„, t)dt+ ^ b ij{ x i,- --, x n, t)dwj 

i=i 

The vector A is called the drift of process X. A is the derivative of the expected value of 
X. The matrix r = BB T is the derivative of the covariance of X. W = (wi, ... ,w n ) is a 
n-dimensional standard Brownian motion. 
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Often, the variables x, are prices of securities available on the market, and are therefore strictly 
positive processes. The expected increments and covariance of increments are then expressed 
in relative values: 

dxi n 
Vi G [l,n], — = (i Xi (xi, ... ,x n ,i)dt+Y^ V(/(*i, ... ,x n ,t)dwj 

Xi 7=1 

with 

[L Xi = di/Xi, Vij = bij/Xi 

The matrix V = { v ij)(ij)e[i,n] 2 * s called volatility matrix, and the covariance of relative returns 
is the matrix 2 = M^n^: 

2 = W T 

We have the following relationship: 

lij — a ij x i x j 

3.2 European securities 

A security is called European security iff future cash flows cannot be influenced by decisions 
from the holder taken after the purchase date (besides of course selling back the security). 
Then the cash flows are only functions of time and information. 

The cash flow generated by a European security C during the time interval dt, assuming C is 
held indefinitely, will be denoted by fc(x(t), t)dt. 

If the price C(t) of a European security C is positive, we can define the instantaneous relative 
cash flow rate or dividend yield as 

d c (x(t),t)=f c {x(t),t)/C{t) 

The aggregate cash flow Fc(to, i) from the purchase date to to any future date t follows the 
equation: 

F c {t 0 , to) = 0, dF c = fcdt = d c Cdt 

The total gain process Gc, sum of the security's liquidative value and the aggregate cash flow 
is: 

G c (to,t) = C(t)+F c (t 0 ,t) 

Let fic denote the expected capital rate of return of C : 

^/dC^ 
Hcdt = E{—) 

The expected total rate of return of C, i.e. the expected rate of return of the total gain process 
is: 

(J-Gc — (J-c+dc 
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The dividend yield of the money market account C, called short term interest rate, is denoted 
by r(x) . If the proceeds of the money market account are continuously reinvested, the total 
gain process L(? 0 , t) follows the equation: 

L(t 0 , to) = 1, dL = r(x)L(to, i)dt 

or equivalently: 

L(t 0 ,t) = exp QT r(x(r))dr^ 

3.3 Arbitrage Pricing 

For any risk factor x,, let e(C/x,) denote the elasticity of price C to x,-: 

(rl \ — Xi dC 

e{c ' Xi) ~ 

The major result of the Arbitrage Pricing Theory is the following. There exist numbers 
Ai , . . . , A„, called market prices for risk, such that for any security C, the following relationship 
holds: 

n 

VGc = r + £<?(C/x,)A,- 

i=l 

Furthermore, if factor x,- is traded, the market price of risk for factor x, is the expected total rate 
of return on x, in excess of the riskless interest rate. In particular, if all x, are traded, we have: 

n n 

VGc ~ r =Y; e ( C / x i)(t 1 G Xi ~r) = ^2e(Clxi)(ii Xi +d Xi - r) 

1=1 !=1 

Using the properties of diffusion processes, the above results lead to the following partial 
differential equation, called Black-Scholes equation: 

~ Tt = {dC ~ r)C + E ^{, Xi - Xfo + - £ ^ 7;7 (2) 

If factor x, is the price of a traded security, the term fi Xj - A, is simply r - d Xj . In particular, if 
all factors are traded, Black-Scholes equation simplifies to: 




We see that the above equation does not depend on the market prices for risk. Therefore, we 
can replace the information process X by the so-called risk-neutral information process X for 
which all market prices for risk are zero. 
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Using a theorem known under the name of Feynman-Kac Formula, we can represent explicitly 
the solution of the above equation: 



where E t represents the expectation under the fictitious risk-neutral information process X. 
Under suitable technical conditions, and using the linearity of the expectation operator: 



In particular, if a security only yields cash flows fc(x(ti), t\), ... ,fc(x(tk), tk) at discrete 
dates t\ < . . . < tk, we have for any time t < t\: 

c{x{l) , t]= p,(mM) 

Therefore, the problem of computing the price of a security with k discrete cash-flows 
fc{x{t\), t\), . . . ,fc{x{tk), tk) reduces to that of computing the prices of each of the k 
securities C h , . . . , C tk with respective single cash flows fc(x(ti), t\ \ . . . ,fc(x{tf c ),tk)- 

In practice, securities only yield cash flows at discrete times. In the following, we will consider 
without loss of generality a security C with a single cash flow fc(x(T)) at a given expiration 
date T. Then, 

«*).0 = *(*gf) 

Let p(x) be the risk-neutral probability density of X(t) knowing X(t). We have: 

c ^=L(0)) p{x)dx (4) 

The practical problem of pricing a European security under n sources of uncertainty reduces 
to that of computing the above n-dimensional integral. 



4 Monte Carlo valuation 

4.1 Numerical integration on the real line 

We consider the problem of evaluating the simple integral: 

1(f) = [ b f{x)dx 

J a 
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on the compact interval [a, b] for any function/. When the above integral cannot be computed 
analytically, the natural solution consists of approximating the integral by a weighted sum of 
values off over a finite number of samples x 1 , . . . ,x m . The simplest quadrature rule is the 
trapezoidal rule: 

;=l 

with 

V?6[l,m-1], x'=a+—(b-a) 

m 

In the above formula, m is an integer chosen large enough for accurate approximation of the 
integral. 

More complex formulas such as Simpson rule or Romberg integration have a higher order of 
accuracy for smooth integrands (see [Stoer and Burlisch, 80] for a review). 

A technique of particular interest is that of Gaussian quadrature. If the function / is the 
product of a polynomial g of degree 2m — 1 and some known kernel W, then we can find points 
x 1 , . . . , x" 1 and some weights w 1 ,.. . ,w m independent of g such that the integral: 

Iw(g)= ff{x)dx= [ b g(x)w(x)dx 

J a J a 

is exactly equal to the sum 

m 

(•=i 

4.2 Multidimensional integration 

The basic idea underlying numerical quadrature generalizes simply to n-dimensional integration 
problems. We consider an n-dimensional hypercube H = [a 1 , b l ] X . . . X [a",b n ] and the 
quadrature problem: 

!{f) = / f{x\ , ■ ■ ■ , x n )dx 1 ...dx n (5) 

For each coordinate x,-, we choose m,- samples This generates a lattice of 

n-dimensional samples containing the mi x ... X m„ points 

VQi,...j„) G [0,mi] x ...[0,m„], X^" = {x 1 ^...,^) 

We can find appropriate weights w /1, --* / " such that the sum: 

]T /(.v^...,v;;y: * 

O'i,-Jn)e[0,mi]x...[0,m„] 

approximates the integral /. For example, if we have for the i th dimension a one-dimensional 
integration formula of the form: 

7=1 
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we can build the Cartesian product rule with weights 

n 

w iu - J " = Y[ 
;-=i 

In particular, if m\ = . . . = m n = m, the number of samples is m". Therefore, evaluating the 
sum is a problem requiring a computation time exponential in the number of dimensions n. 
Hence, direct numerical integration is impractical for high-dimensional problems. 

Practical computation of a high-dimensional integral requires to choose fewer samples 
X 1 , . . . , X M and weights w 1 , . . . , v/ 4 , such that the total number of samples M is polynomial in 
the dimension n. If the samples and weights are chosen carefully, the sum 

M 

!=1 

approximates the integral /. 

Choosing appropriate samples and weights is in general a difficult (i.e. intractable) problem. 
However, in many practical situations, the functions/ to be integrated have specific properties 
that make the problem tractable. In particular, security pricing problems of the type (4) lend 
themselves surprisingly well to numerical integration. 

The main idea described above underlying Gaussian quadrature can be extended to multi- 
dimensional spaces. While classical Cartesian product rules require a number of samples 
exponential in the dimension, there exist non-Cartesian rules requiring only a polynomial 
number of samples. Let us consider a kernel W(xi ,x n ) defined over R". In general, it can 
be shown (see [Stroud, 71]) that there exist a Gaussian rule exact for all functions/ that can 
be represented as the product of W with a polynomial of degree m, and that requires at most 



samples, number which is indeed polynomial in n for a given m. For example, there exist 
a formula exact for polynomials of degree two in R" requiring only M = (n + l)(n + 2)/2 
samples. In fact, for second degree polynomials, there even exist a formula with only n + 1 
samples in dimension n. In practice however, those rules have some negative coefficients w', 
and are thus very unstable numerically. Most quadrature rules with only positive coefficients 
require at least 2" samples, and thus cannot be used for very high dimensional spaces. 
[Stroud, 71] presents a comprehensive review of results on generalized Gaussian quadrature 
rules. 

We experimented with the best rules advocated by [Stroud, 7 1 ] on typical asset pricing problems 
in dimensions ranging from 1 to 10. In particular, we used the Cartesian product Gauss rule, 
exact for any polynomial of degree 3 multiplied by a Gaussian kernel W(X) = exp(- j||X|| 2 ). 
This rule requires 3" samples. In dimension 10, this represents almost 60,000 samples. 
Although the rule is exact for any polynomial of degree 3, our experiments have shown (see 
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Section 6) that it performs very poorly of typical asset pricing problems. We also used the 
spherical Lobatto rule, exact for any polynomial of degree 5, and requiring 2" +1 — 1 = 2047 
samples in dimension 10. The results, presented in Section 6, show that these Gaussian 
quadrature rules are inadequate for asset pricing problems in more than a few dimensions. 

A alternative way of choosing the samples is to generate them randomly, using a pseudo- 
random number generator. The resulting integration technique, known under the name of 
Monte Carlo integration, is quite widespread and simple to implement. We present in the next 
subsection the general principle underlying Monte Carlo integration. 

4.3 Monte Carlo method and importance sampling 

We consider the hypercube H = [0, 1]" without loss of generality, and a uniformly distributed 
random vector X over H. By definition, the expected value of the variable/(x) is the integral 
/(/) defined by equation (5): 

E U (/(X)) = 1(f) = [ f( Xl x n )dx x ...dx n 

where E u denotes the expectation under the uniform distribution. The principle of Monte Carlo 
integration is to generate M independent random samples X 1 , . . . , X M of the vector X and to 
compute the empirical mean: 

1=1 

By the law of large numbers, this sample mean is approximately equal to the integral /(/) for 
M large enough. 

This above crude Monte Carlo method converges very slowly even in dimension 1 , and usually 
does not converge in a reasonable amount of time in high-dimensional spaces. 

One technique for accelerating convergence consists of performing an appropriate change of 
variables in the integral. The technique is known in the literature as importance sampling. 
Let us assume that we know a positive function/ over H, such that/ approximates/ in some 
sense, and that the integral: 

= / , • • • , x n )dxi ...dx n 
Jh 

is easy to compute. Then, we can consider the probability distribution whose density is 

_ / 
q ~W) 

Instead of generating uniformly distributed samples, we generate the samples 

X 1 , . . . , X M with the importance sampling density q. To do so, we assume that we can find a 

change of variable in R" Y = (vi , . . . , y„) such that: 

\det(dY/dX)\ = q 
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Therefore, 

/(/) = ffdX= [ f -qdX = [ f -dY = &{ f -) 
Jh Jh q Jy(h) q q 

where E q denotes the expectation for the distribution of density q. Equivalently: 

/(/) = KfMp 

In particular, when / is positive, if we choose / = /, /// is always equal to 1, hence 
£?(///) = 1 and there is no need to generate any samples. However, the normalizing factor 
l(f) = /(/) is precisely what we aimed at estimating in the first place. In practice, / must be 
chosen in such a way that it is approximately proportional to/, while still significantly easier 
to integrate than /. 

There is no easy systematic way of choosing the appropriate importance sampling density 
q. However, the integration problem is often naturally expressed in such a way that / is 
proportional to a special probability density p: 

f = 8P 

For problem (5), p is the function constant over the hypercube H, density of the uniform 
distribution. In general, when the problem at hand is to compute the expected value of a 
random variable g(x) where X has density p, we have / = gp. p is then called the natural 
importance sampling density, and the method of importance sampling simply consists in 
simulating X with its natural probability distribution p. 

For asset pricing problems of the form (4), the risk-neutral probability density p{X) of the 
information process X is often a good candidate for importance sampling. We generate samples 
X 1 , . . . , X M with the risk neutral probability distribution, and approximate the asset price by 
the average discounted cash flow: 

This technique is called risk-neutral importance sampling. 

Recently, [Dantzig and Glynn, 89] and [Dantzig and Inf anger, 91] have proposed the method 
of additive importance sampling, an interesting paradigm for choosing importance sampling 
functions. They consider a natural probability density function p(x) in R n , which makes the n 
coordinates independent two by two: 

Pi*) = f[Pi(*i) 
!=i 

The corresponding distribution could be for example the uniform distribution over an hypercube 
H(Pi = l/(ft,-a,))orthemultivariatestandard normal distribution (p,- = exp(-x?/2)/\/27r)- 
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We consider the problem of computing the integral: 

h(f) = I f{x)p{x)dx 

The method of additive importance sampling consists of identifying n positive functions 
/i(*i), . . . ,f„(x„) such that their sum 

n 
i=l 

approximate the function/. The integral off is easier to compute than that off, since it is a 
sum of monodimensional integrals: 

" P+OO 

W) = E / fix)pi{x)dx 

i=l J -°° 

Then 

Ip ( / ) = J/f P , x =±J/M fl{xi)p{x)dx 

The problem reduces to that of computing each of the n integrals: 



V/ 6 [1,72], /' = 




Due to the particular structure of the above integral, it is easy to generate samples with a 
density proportional to/(x,)p(x). Then, iff is actually a good approximation off, the additive 
importance sampler will converge substantially faster than the simple natural importance 
sampler. The importance sampling density is 

q=pf/I P (f) 

It is the product of an additive function and of a multiplicative function of the variables 
X\ , . . . , x n . 

4.4 Quasi-random sampling 

Monte Carlo methods combined with appropriate importance sampling schemes are exper- 
imentally efficient. However, generating the samples randomly is not necessarily the best 
way of obtaining an accurate quadrature rule. Indeed, the efficiency of Monte Carlo methods 
does not depend on the random nature of the samples, but on their equidistribution properties. 
Deterministic techniques for generating samples often perform better than random techniques. 
For historical reasons, deterministic techniques for choosing the samples are called Quasi 
Monte Carlo techniques. There is a fairly large literature on Quasi-Monte Carlo methods, that 
we will not present here. We refer the reader to [Zaremba, 68] and [Haber, 70] for a review. 
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We developed a very simple method for generating deterministically equidistributed samples, 
that happens to perform very well in practice. For low dimensional problems (say, problems 
in dimensions < 5), this approach is considerably more efficient than classical methods based 
upon random number generation. On higher dimensional problems, the method yields results 
similar to that of classical Monte Carlo methods. Hence, we can get the best of both worlds: 
reasonable efficiency for high-dimensional problems, and high efficiency for low dimensional 
problems. 

We consider again problem (5). We assume that an appropriate importance sampling density q 
has been chosen. After the corresponding change of variable Y(X) such that \det(dY/dX)\ = q, 
we are faced with the following problem: 

l{g) = [ gdY 

with g = f I q. Instead of generating randomly uniformly distributed samples, we choose 
to sample the ^-dimensional hypercube deterministically. We first choose an appropriate 
quantization step h for each coordinate, say h = 1 /2 k for a positive integer k. Since we make 
our calculations on a 32 bits computer, a simple choice is k = 32. A quantized point Y in the 
hypercube is then represented by a sequence of n £-bits integers (yi , . . . , y„). We consider the 
Toi-bits integer: 

L(y u ...,y n ) = £ yi 2 k ^) 
i=i 

For each number y,-, we can write the representation of y ( - in base 2: 

» = 1 

7=1 

where y| G {0, 1}. Then: 

L(y u ...,y n ) = ±j:^-^ 
i=i j=i 

The mapping L is a one-to one mapping from the quantized hypercube [0, 2 k [" onto the integer 
interval [0, 2 nk [. Intuitively, this mapping corresponds to scanning the quantized hypercube in 
lexicographic order. 

Alternatively, we can consider the mapping in Peano order, defined as: 

4 .) = EE^- |)+M 

7=1 (=1 

If we consider the variables 

V/6[l,*], : i = Y J ^ ' 
i=i 
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spanning the hypercube [0, 2 n [ k , 



we get: 



k 



P(yu...,yn) = Y,zA n{j - 1] 



Hence, the Peano mapping can be viewed geometrically as a one-to-one mapping between the 
hypercube [0, 2 [" and the hypercube [0, 2"[ . In our experiments, the lexicographic mapping 
L and the Peano mapping P gave similar results. For n = 1, both the lexicographic and the 
Peano scan reduce to the identity function over the interval [0, 2 k [. 

Given any of those two mappings, say L, we can now construct a sampling of the hypercube 
by first sampling regularly the interval [0, 2 nk [, and then taking the inverse of the mapping. 

Given the desired number of samples M, we perform the integer division: 

2 nk = MA + B 

with B < M. Then, we sample regularly the interval [0, 2 nk [ with the M numbers 

Vr£[l,M], Y r = rA 

The corresponding samples in the hypercube [0, 2 k [ n are the M points (y[ , . . . , y r n ) verifying 



If M is properly chosen, the samples are correctly distributed over the hypercube. In practice, 
we choose M of the form M = p a , where p is a prime integer different from 2, say p = 3. 
More advanced number theoretic methods of this type are described in [Stroud, 71]. 

In dimension n = 1, the sampling is simply a regular sampling of the interval [0, 1]. 
Hence, the quadrature formula is almost identical to a regular trapezoidal formula, much 
more accurate than a formula generated by a random sampling. More generally, this quasi 
Monte Carlo method has the advantage of not introducing spurious randomness in low 
dimensional integration problems. For high-dimensional problems, experiments have shown 
that it performs similarly to classical Monte Carlo methods. 

5 Quadratic resampling 

We consider again the problem of computing the integral: 



where p is a natural probability density function, and/ any integrable function. In particular, 
the general quadrature problem (5) is a special case of the above problem, p being the uniform 
density over the hypercube H. If we consider the random vector X with density p, this is 
equivalent to the problem of computing the expected value 



L{y\,...,y r n ) = Y r = rA 



W) = [ f{x)p(x)dx 

JR" 



I p (f) = E(f{x)) 
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We assume that the n first order moments have been precomputed (either analytically or 
through Monte Carlo simulation, using one of the importance sampling techniques described 
above). 

V/ 6 [1,72], l ( Xi )=E(x i )= [ x iP (x)dX 

The vector whose coordinates are £(x,) is denoted E(X), and called expected value of X. We 
also assume that the n 2 second order moments have been precomputed. 

V(jj') G [l,n] 2 , I p {xiXj) = E{xiXj) = / XiXjp{x)dX 

jRn 

The matrix of second order moments is denoted E(XX T ). By definition, the covariance matrix 
ofXis: 

k x = e((x- e(x))(x - e(x)) t ) = e(xx t ) - e(x)e(x) t 

Kx is a non-negative symmetric matrix. Without loss of generality, we will assume that K is 
positive. Indeed, if K x is degenerate, we can always diagonalize it in an orthonormal basis, 
perform the corresponding change of variables in the integral I p , and discard all degenerate 
dimensions. Then, the above quadrature problem is transformed in a lower dimensional 
non-degenerate quadrature problem. 



Since K x is non-negative symmetric, we can find a square root, i.e. a matrix denoted <jKx 
such that: 

I — I — T 
K x = yKxyKx 

This can be done either through LU decomposition, or through diagonalization of Kx in a 
orthonormal basis. Since K x is positive, \fK~x is regular (i.e. has an inverse). 

We then consider a sequences of samples X 1 , . . . , X M , generated either deterministically or 
randomly. For example, they can be generated randomly so as to follow the distribution 
of density p. In general, we can choose any appropriate importance sampling density q for 
generating the samples. Then, we can define the empirical mean approximating the integral 
I p (f) as: 

M 



f(x) = £/(*V 
(•=1 

with 

V/ G [l,n], W l = 77 7 T 7 T" —, 2 

EjLip{^)/q{Xj) q{x) 
The normalizing factor for w"s has been chosen so as to have: 

M 



In particular, the empirical mean is: 



Em/ 

!=1 



M 

X = Y J X i w i 

!=1 
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We can define similarly the empirical covariance: 



K X =(X- X)(X - X) T = XX T -XX T 



More explicitly: 



M M 

= - x){x' - X) V = ^Z'X'V -If 
(=1 (=1 

If the samples are generated randomly, the law of large numbers applies. Therefore, for a 
sufficiently large number M of samples, the empirical mean X and the empirical covariance 
Kx will be arbitrarily close to the real mean E(x) and the real covariance Kx- Therefore, 
since the set of regular matrices is open, the empirical covariance Kx must be non-singular 
for M sufficiently large. If the samples are not generated randomly, there exist deterministic 
sampling techniques guaranteeing the same convergence property of the empirical mean and 
covariance to the real mean and covariance, and the same result applies. 

Since Kx is non-singular, its square root \fKx is also non-singular. 

We define the gain matrix 



-l 

H = \/Kx\l Kx 



and the new random variable: 



Y = H(X -X)+ E(X) 
We consider the M samples Y l = H(X - X) + E(X). For this particular sampling, we get: 

M 

Y=J2 Y'w' = E(X) 
i=i 

Similarly: 

K^= (Y -Y)(Y -Y) t = (H(X -X))(h(X -X)) t = HKx~H T 
Using the expressions above for H and Kx as products of square root matrices we get: 

Ky = Kx 

Hence, the empirical first and second order moments using the samples Y 1 are exactly equal 
to the real first and second order moments of X. In particular, the empirical mean of any 
polynomial/ of degree two or less in the variables x\ , . . . , x„ verifies: 

M 



f(Y) = J2f(Y i )w i = E(f(x)) 



We can state: 
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Any numerical quadrature formula generated through quadratic resampling is 
exact for any polynomial of degree two or less 

In the above sense, the method of quadratic resampling can be viewed as a particular Gaussian 
integration rule of degree 2, following the terminology of [Stroud, 71]. 

The method of quadratic resampling can be summarized as follows: 

1) Precompute the expected value E(x) and covariance matrix Kx of the random vector X 

with density p, if they are not available in analytical form. 

2) Choose an importance sampling density q(x) well-suited to the problem at hand, the 

simplest choice being q = p. 

3) Sample the random vector X using any available technique, i.e. Monte Carlo sampling 

or deterministic sampling, to obtain a sequence X 1 , . . . , X M . Define the empirical 
quadrature formula: 

M 

f(x) = YA*W 

Compute the empirical mean X and empirical covariance matrix Kx. Choose M large 
enough so that Kx is non-singular. 

4) Compute the matrix H = \fKx\fKx 

5) Compute the modified samples T = H(X - X) + E(X) , to obtain the modified empirical 

quadrature formula: 

M 

f(Y) = £/(r'V 

!=1 

6 Experimental results 
6.1 A test case 

We implemented the different numerical quadrature schemes presented above on a test case 
described below. 

We consider two real numbers a and p with — l/(n — l) < p < 1, and the non-negative 
symmetric matrix^ = (^>)(, i/ ) e [i,„]2, such that 

Vi G [l,n], hi = o- 2 

and 

V(ij') G [l,n] 2 , i 7~ ji kj = P° 2 
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We consider a random vector X = (xi,...,x n ) T G R n , following a multivariate normal 
distribution with zero mean and covariance matrix K. 

Let us assume that K is of rank k < n, and let us compute a square root V of K, i.e. an x k 
matrix such that: 

K=W T 

We consider the random vector Z: 

Z = (v T v)~ 1 v T x 

By construction, X = VZ, and Z follows a k-dimensional standard normal distribution. 

We can always take k = n by completing the matrix V with columns of zeros, and consider 
the quadrature problem: 

l(f) = E(f(x))= [ f{x) Px {x)dX= f f{VZ)p z {z)dZ 

where px and pz are the densities of X and Z. Hence: 



Kf) = J /(t ^*pW/2) JZ 

•/ R" V27T 



We choose the particular expression for/: 

/(xi, . . . ,x n ) = maxexp(x,) 
(=i 

With this choice, we obtain a three parameters quadrature problem: 

I{n,a,p)= / max exp( > v^J ^ „ — ^<fei . . . dz„ 

JR" '=1 ~] V27T 

The above problem admits the following financial interpretation. We consider a stock market 
including n dividend-free stocks with prices S\,...,S n following a jointly lognormal diffusion 
process with any arbitrary expected yearly returns and a covariance matrix of yearly returns 
equal to K. We have: 

Vi G [1 , n], = mdt + v v dw j 

* 7=1 

where W = {yv\, ... , w n ) T is an n-dimensional standard Brownian motion. It can be shown 
that the solution of the above equation is 

n 

V/ G [l,n], Si{t) = Si{0) exp((/x,- - k u /2)t + £ v ijWj (t)) 

7=1 

Using ku = a 2 and setting x,(f) = J2j=i v ij w ji. t \ we § et: 

Si(t) = Si(0) exp((/x,- - <t 2 /2» exp(x,(0) 
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We assume that the annualized interest rate r is constant and known at time 0. Therefore, 
the risk-neutral information process is obtained by replacing all /x, 's in the above equations 
by r. If we consider the price C at time zero of a contingent security with a single cash flow 
g(Si, ... ,S„) at date T, the Feynman-Kac formula yields: 

C= [ exp(-iT)g( , . . . , S n ) 6XP( " gfi Z * /2) dzi... dz n (6) 

JR" V27T 

with 

n 

Si = 5,(0) exp((r - ct 2 /2)t) exp((^ v ijZj ) Vf) 

7=1 

If the contingent claim is an option with zero strike entitling its holder to choose any one of 
the n securities S\, . . . , S n at expiration date, we have: 

g(S 1 ,...,S n ) = max Si 
i=i 

Indeed, any clever investor will choose the security with the highest market value at time T. 

Finally, if we assume that the initial prices of the stocks at time zero are all equal to 1 , and that 
the cash flow is exactly one year hence (T = 1.0), we get: 

Si = exp(r - <r 2 /2) exp(x,) 

Hence: 

C = exp(-<7 2 /2)/(?i, a, p) 

The above quadrature problem l[n, a, p) can be interpreted as an option pricing problem. 
Although this problem can be solved analytically using the ^-dimensional multivariate normal 
distribution function (see [Johnson, 87]), it is an appropriate testbed for numerical integration 
methods. 

We tested four different Monte Carlo procedures, to be described below, on the above test 
case for different values of the parameters. All the four methods use importance sampling. 
Indeed, for large values of n, a crude Monte Carlo method using uniformly distributed 
samples z.\ , . . . , z n would take an astronomical number of samples to converge. The simplest 
importance sampler is the Natural Importance Sampler (NIS), which consists of generating 
Z\ , . . . , z n following a n-dimensional standard normal distribution. Since the normal distribution 
has the property of central symmetry at the origin, we also use in the four methods the technique 
of Antithetic Variables (AV). Instead of generating randomly M samples, we generate M/2 
samples Z 1 , . . . , z M / 2 and take for the remaining samples z M / 2+1 , . . . , Z M the opposite values 
— Z 1 , . . . , — Z M / 2 (M is assumed even for simplicity). The technique of antithetic variables 
applied to normal samples guarantees in particular that the empirical first order moments (i.e. 
the expected value of Z) are exactly equal to their theoretical values (0 here). If the function/ 
has some amount of correlation with the coordinates zu the technique will reduce the variance 
of the samples, and thus improve the accuracy of the Monte Carlo estimate. This is the case 
in particular for the function/ chosen in our test case. In general, it is the case for most asset 
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pricing problems, since the cash-flows generated by securities are often monotonous functions 
of the underlying economic factors. 

We now describe the four methods: 

Method 1: AV-NIS Method 1 is simply the one just described. It uses the technique of 
antithetic variables, and uses the natural importance sampler. 

Method 2: AV-NIS-AIS Method 2 uses in addition to Method 1 the Additive Importance 
Sampler (AIS) described in the previous sections. The functions/} chosen are: 

V/ G [1,72], fi(zi) = exp(v/iZ/) 

Method 3: AV-NIS-QR Method 3 uses in addition to Method 1 the technique of Quadratic 
Resampling (QR). Since the moments of the multivariate normal distribution are known 
analytically, we do not need to precompute E(Z) and Kz'. 

E(Z) =0, K z = I n 

where /„ is the identity matrix in R n . 

Method 4: AV-NIS-AIS-QR Method 4 uses in addition to Method 2 the technique of 
Quadratic Resampling. 

These four methods were tested on the quadrature problem l(n, a, p) for different values of 
the parameters. All tests used a sample size M = 4000. Table 1 shows the standard errors for 
the four methods in dimension n = 10 using a = 1 and four different values of the correlation 
parameter p = —1/9,0, 1/2, 1. The value a = 1 corresponds to a volatility of 100 percent 
per year for each of the 10 prices Si,... , Sio- 

Table 2 shows the standards errors, again in dimensions = 10, using a = 0.1. This represents 
a volatility of 10 percent per year for each of the underlying assets prices. 

The results presented in Tables 1 and 2 illustrate the relative efficiency of the different 
methods on 10-dimensional problems. For problems with high volatility (<r = 100%), 
quadratic resampling (Method 3, AV-NIS-QR) reduces the error by a factor ranging from 2 
(low correlation) to 3 (high correlation) (as compared to Method 1, AV-NIS). The method of 
additive importance sampling also reduces the error by a factor of 2 (low correlation) to 1 .2 
(high correlation) (Method 2 / Method 1). Interestingly, the error reduction power of quadratic 
resampling and additive importance sampling can be combined (in Method 4), yielding an 
error reduction factor of approximately 4 over Method 1. This corresponds in terms of 
computing time to a speedup of 16 for a given accuracy. The quadratic resampling method 
works best on problems exhibiting a high correlation between variables p > 0.5, while the 
additive importance sampling method works best on problems with low correlation between 
variables. 
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p 


7(10,1.0^) 


Method 1 


Method 2 


Method 3 


Method 4 


-1/9 


5.91 


0.0734 


0.0421 


0.0423 


0.0260 


0.0 


5.62 


0.0697 


0.0311 


0.0424 


0.0204 


0.5 


4.18 


0.0534 


0.0219 


0.0287 


0.0125 


1.0 


1.65 


0.0265 


0.0213 


0.0083 


0.0054 



Table 1: Standard errors in dimension n = 10 for a = 100% 



p 


/(I0,0.1,p) 


Method 1 


Method 2 


Method 3 


Method 4 


-1/9 


1.178 


0.00109 


0.00104 


0.00050 


0.00051 


0.0 


1.168 


0.00102 


0.00097 


0.00047 


0.00047 


0.5 


1.119 


0.00070 


0.00065 


0.00032 


0.00032 


1.0 


1.005013 


0.00016 


0.00013 


0.00000048 


0.00000032 



Table 2: Standard errors in dimension n = 10 for a = 10% 

For lower volatilities (<r = 10%), the additive importance sampler performs poorly. The 
method of quadratic resampling yields an error reduction factor ranging from 2 for problems 
with low correlation to several hundreds for problems with high correlation. Indeed, the 
cash-flow function is closer to its second order Taylor expansion for lower volatilities. Since 
quadratic resampling yields exact results for functions that are equal to their second order 
Taylor expansion, it is likely to give more accurate results for problems with low volatilities. 

All calculations were performed on a DECstation 5000 model 200 desktop workstation, using 
a MIPS R3000 microprocessor at a clock rate of 25 Mhz. All calculations were performed 
in double precision floating point arithmetics. The simulation program was written in the 
C programming language, and compiled without any optimization. The resulting computing 
times for each run of 4000 samples were just below 1 second. 

Tables 3 and 4 present the same results in dimension n = 5. The conclusions are similar. 
For high volatilities, quadratic resampling reduces the error by a factor ranging from 2 (low 
correlation) to 5 (high correlation). When combined with additive importance sampling, 
the total error reduction factor ranges from 5 to 10, corresponding to a reduction factor in 
computing time ranging from 25 to 100. For lower volatilities, the additive importance sampler 
performs again poorly. The method of quadratic resampling yields an error reduction factor 
ranging from 5 for problems with low correlation to several hundreds for problems with high 
correlation. In Table 4, the error reduction of Method 4 over Method 1 ranges from 5 to 1000, 
corresponding to a reduction factor in computing time ranging from 25 to 1 million. The 
resulting computing times for each run of 4000 samples were about 1/3 of a second. 

Finally, Tables 5 and 6 present results in dimension n = 100. For high volatilities (Table 
5), quadratic resampling reduces the error by a factor ranging from 1.25 (low correlation) to 
3.5 (high correlation). When combined with additive importance sampling, the total error 
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p 


1(5,1.0, p) 


Method 1 


Method 2 


Method 3 


Method 4 


-1/4 


4.39 


0.0719 


0.0365 


0.0230 


0.0141 


0.0 


4.07 


0.0594 


0.0186 


0.0255 


0.0088 


0.5 


3.29 


0.0448 


0.0102 


0.0205 


0.0058 


1.0 


1.65 


0.0284 


0.0183 


0.0067 


0.0028 



Table 3: Standard errors for n = 5 and a = 100% 



p 


I(5,0A, p) 


Method 1 


Method 2 


Method 3 


Method 4 


-1/4 


1.141 


0.00133 


0.00126 


0.00027 


0.00027 


0.0 


1.126 


0.00117 


0.00110 


0.00024 


0.00024 


0.5 


1.090 


0.00079 


0.00073 


0.00017 


0.00016 


1.0 


1.005013 


0.00017 


0.00011 


0.00000040 


0.00000017 



Table 4: Standard errors for n = 5 and a = 10% 



p 


7(100, 1.0,/?) 


Method 1 


Method 2 


Method 3 


Method 4 


0.0 


13.58 


0.1235 


0.0984 


0.0994 


0.0818 


0.5 


7.94 


0.0850 


0.0694 


0.0486 


0.0422 


1.0 


1.65 


0.0301 


0.0294 


0.0087 


0.0084 



Table 5: Standard errors for n = 100 and a = 100% 



p 


7(100, 0.1, p) 


Method 1 


Method 2 


Method 3 


Method 4 


0.0 


1.286 


0.000901 


0.000887 


0.000698 


0.000696 


0.5 


1.1975 


0.000588 


0.000579 


0.000457 


0.000456 


1.0 


1.005013 


0.000169 


0.000165 


0.00000051 


0.00000049 



Table 6: Standard errors for n = 100 and a = 10% 
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reduction factor ranges from 1.5 to 3.6, corresponding to a reduction factor in computing time 
ranging from 2.3 to 13. For lower volatilities, the additive importance sampler performs again 
poorly. The method of quadratic resampling yields an error reduction factor ranging from 1.3 
for problems with low correlation to several hundreds for problems with high correlation. In 
Table 6, the error reduction of Method 4 over Method 1 ranges from 1.3 to 350, corresponding 
to a reduction factor in computing time ranging from 1.7 to 100,000. 

The resulting computing times for each run of 4000 samples were about 1 minute and 
10 seconds. In general, the computation time is quadratic in the dimension n, since the 
computation of each sample requires to compute the transformation X = VZ, and V is a n x n 
matrix. 

To illustrate typical computing time in practical cases, let us consider the computation of 
7(10,1.0,0.0) using Method 4. The result, presented in Table 1, is 5.62, and has a standard 
error of approximately 0.02. Using the central limit theorem, we see that the accuracy for 
a confidence level of 99.995% is approximately 0.08 (4 times the standard error). In order 
to make the last digit significant, an accuracy of 0.01 (8 times higher) is required. Hence, 
the computation time for 0.01 accuracy would be 8 X 8 = 64 times higher. The result 
would be obtained in about 1 minute of computation on the same hardware platform, with 
256000 samples. With Method 1, the same result would be obtained in about 12 minutes of 
computation, using more than 2.5 millions of samples. Not surprisingly, Table 2 shows that 
problems with lower volatility require much less computing time. Even for low correlations, 
three digit accuracy is achieved in a small fraction of a second, and four digit accuracy in a 
couple of seconds (the error margin of Method 3 for p = 0.0 is about 4 x 0.00047 & 0.002 for 
a confidence level of 99.995% using only 4000 samples, i.e. a computation time of 1 second). 

In general, the computing time may vary from several minutes to a fraction of a second, 
depending of the characteristics of the covariance matrix, and on the accuracy required. Since 
the accuracy of Monte Carlo methods is proportional to the square root of the number of 
samples, getting one additional significant digit requires a computing time 100 times longer 
on a given hardware platform. 

6.2 Comparison with classical Gaussian rules 

We implemented several classical Gaussian integration rules proposed in [Stroud, 71] and 
experimented with them on the above test case. 

We present below some results for 2 different rules, designed for approximating integrals of 
the form: 



Gaussl is a classical Cartesian product rule of the form described in Section 4.2. It is exact if 
/ is a polynomial of degree 3. 

In dimension 1, the three sample points are -\/3, 0, and \/3. The corresponding weights 
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p 


7(10,1.0^) 


Gauss 1 


Gauss2 


Stdev Method 4 


-1/9 


5.91 


5.53 


5.46 


0.034 


0.0 


5.62 


4.901 


5.35 


0.028 


0.5 


4.18 


4.11 


3.64 


0.018 


1.0 


1.65 


1.64 


1.64 


0.0078 



Table 7: Values computed from Gaussian quadrature rules (n = 10, a = 100%) 



p 


7(10,0.1^) 


Gauss 1 


Gauss2 


Stdev Method 4 


-1/9 


1.178 


1.180 


1.165 


0.00065 


0.0 


1.168 


1.159 


1.142 


0.00061 


0.5 


1.119 


1.120 


1.096 


0.00041 


1.0 


1.005013 


1.005013 


1.005013 


0.0000004 



Table 8: Values computed from Gaussian quadrature rules (n = 10, a = 10%) 

are 1 /6, 2/3, and 1 /6. In dimension n, the total number of samples is therefore 3". 

Gauss2 is a so-called spherical Lobatto rule (see [Stroud, 71]). It is exact for polynomials of 
degree 5. The total number of samples in dimension n is 2" +1 — 1 

In dimension 10, the number of samples required for the Cartesian product rule Gauss 1 is 
59049. The number of samples required for the spherical Lobatto rule Gauss2 is 2047. Results 
are presented in Tables 7 and 8. In order to compare these Gaussian rules with the Monte 
Carlo method, the last column of the tables shows the standard error of Method 4 above, for 
the exact same number of samples as the Lobatto rule, i.e. 2047 in this case. The results 
from Table 7 show that both Gaussian rules are very inaccurate for pricing problems with high 
volatilities and low correlations. The error is in most cases over 5%, and exceeds 10% on 
several examples, while the standard error of Method 4 remains under 0.5% for 2047 samples. 
Correct results are obtained only on problems with very high correlation. Results in Table 
8 indicate that Gaussian rules perform better on problems with low volatilities, but are still 
not nearly as good as Monte Carlo methods for a given number of samples. The error of the 
Lobatto rule ranges from 1 % to over 2%, while the error of Method 4 remains below 0.05% 
with the exact same number of samples. The Cartesian product Gauss rule Gauss 1 is about as 
accurate as Method 4 for low volatilities, but it requires 30 times more samples. 

Tables 9 and 10 present results in dimension n = 5. Gauss 1 requires 243 samples, while 
Gauss2 requires 63 samples. 

Results are presented in Tables 9 and 10. As before, the last column of the tables shows 
the standard error of Method 4 above, for the exact same number of samples as the Lobatto 
rule, i.e. 63 in this case. The error reported in Table 9 (high volatilities) ranges from 5% 
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p 


1(5,1.0, p) 


Gauss 1 


Gauss2 


Stdev Method 4 


-1/4 


4.39 


4.30 


4.35 


0.116 


0.0 


4.07 


3.78 


3.87 


0.077 


0.5 


3.29 


3.30 


2.93 


0.048 


1.0 


1.65 


1.64 


1.64 


0.020 



Table 9: Values computed from Gaussian quadrature rules (n = 5, a = 100%) 



p 


I(5,0A, p) 


Gauss 1 


Gauss2 


Stdev Method 4 


-1/4 


1.141 


1.115 


1.111 


0.0025 


0.0 


1.126 


1.113 


1.108 


0.0022 


0.5 


1.090 


1.086 


1.073 


0.0015 


1.0 


1.005013 


1.005013 


1.005013 


0.0000001 



Table 10: Values computed from Gaussian quadrature rules (n = 5, a = 10%) 

to 10% on most examples, while the standard error of Method 4 remains under 2% for 63 
samples. Results in Table 10 indicate that Gaussian rules perform better on problems with 
low volatilities, but are still not nearly as good as Monte Carlo methods for a given number of 
samples. The error of the Lobatto rule ranges from 1% to 2%, while the error of Method 4 
remains below 0.22% with the exact same number of samples. 

Of course, these Gaussian rules cannot be implemented in dimension n = 100, since they 
would require an astronomical number (over 10 30 ) of samples. 

In conclusion, we see that classical Gaussian quadrature rules do not stand the comparison 
with the best Monte Carlo methods, even on relatively low-dimensional problems (e.g. n = 5). 
Gaussian integration methods are only useful on problems with very few dimensions (say 1 to 
3). 

6.3 Massively parallel implementation 

We implemented the Monte Carlo procedure on a DECmpp 12000 massively parallel deskside 
workstation with 4K (4096) processors. We present below comparative computing times with 
a DECstation 5000 model 200 workstation. 

We consider again n stocks with a covariance matrix of relative returns of the form described 
in the previous subsection. The price of any claim contingent to the n stocks with a single 
cash-flow date is computed from formula (6). 

We considered several examples for the cash-flow function g. The first example was that of 
an option on the best performance of the n stocks, with a non-zero strike price H. Then, the 
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function g is: 

g(Si , . . . , S„) = max (o, max 5, - h) (7) 

The second example was that of an option on a market index (CAC 40) representing the 
weighted average of n = 40 stocks: 

(40 
0, 
i=i 

Several other forms of cash-flow functions were also tested such as the quality delivery option 
of futures contracts: 

g(S u ...,S n ) = minS ; 

i=i 

We implemented for the above pricing problem the Method 1, AV-NIS, on the DECmpp. 
The program was written in the MASPAR MPL language, a parallel extension of C. The 
parallelization paradigm for Monte Carlo quadrature algorithms is particularly simple on such 
SIMD (Single Instruction Multiple Data) computers. It consists of distributing on different 
processors the computation of different samples, and finally adding the results obtained by all 
processors. We give below experimental computing times for the first example of cash flow 
function. The parallelization speedups obtained for all other examples were strictly similar. 
In this first example using function (7), the initial prices of the n = 10 stocks were all set to 
5 f (0) = 40.0 dollars, the strike price was H = $40.0, the interest rate was set to r = 7%, 
the time to expiration T = 200 days. The volatility a was set to 30%, and the correlation p 
to 50%. In order to reach a 1/2 cent accuracy at a confidence level of 99.995%, we used 16 
millions of samples. The computed price was C = $12.26. 

The computation time on the 4K DECmpp was 62 seconds, using single precision floating 
point arithmetics. In comparison, the computation time for the same problem on a DECstation 
5000-200 in single precision and using all available compiler optimizations was 31 minutes. 
This represents a speedup factor of 30. In double precision, the speedup was 23. We estimate 
that the speedup on a full-size (16K) DECmpp 12000 would be exactly four times the speedup 
obtained on the 4K DECmpp, since the Monte Carlo method can be fully parallelized. Then, 
the speedup factor would be about 90 in double precision and 120 in single precision. These 
figures are clearly not specific to this particular example, and apply in fact to any Monte 
Carlo integration problem where the number of samples required is higher than the number of 
available processors. 

We can conclude that the massively parallel implementation yields a speedup factor of roughly 
two orders of magnitude for most Monte Carlo integration problems. 

Method 3, AV-NIS-QR, parallelizes exactly in the same fashion after precomputation of the 
modified samples. For the particular example above, the error reduction factor using quadratic 
resampling was 2.3, yielding a speedup factor of more than 5 over Method 1, AV-NIS. 

We also compared Method 1 and Method 3 on the same example in dimension n = 100, all 
other parameters being unchanged. The call price is $20.23. The speedup factor obtained 
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through quadratic resampling is only 1.5 in that case. Using Method 3 with 4000 samples, the 
standard error is 0.053, and the computation time is about 1 minute on a DECstation 5000. 
The same computation on a 4K DECmpp 12000 would be done in about 2 seconds. 

Finally, we compared Method 1 and Method 3 on the same example in dimension n = 5. The 
call price is C = 9.70. The speedup factor obtained through quadratic resampling is 12. Using 
Method 3 with 4000 samples, the standard error is 0.024. The computation time is about 1/3 
of a second on the DECstation, and could be done in 1/100 of a second on the DECmpp. 

7 Conclusion 

In this article, we described a systematic numerical technique for pricing arbitrarily complex 
European securities. Besides its obvious applications to trading and hedging in organized 
capital markets, European security pricing has many important applications in various areas 
of risk management such as assets and liabilities management and corporate investment 
decision making. Using this technique, we were able to compute the prices of the most 
complex European instruments in a few tens of seconds on a workstation, and within a 
fraction of a second on a massively parallel computer. We have shown that this technique 
outperforms previous general purpose European security pricing methods, either in the range 
of applicability, or in the accuracy obtained for a given computation time. 

Our approach essentially relies on the well known Monte Carlo method, combined with 
appropriate importance sampling schemes and a new error reduction technique, that we call 
quadratic resampling. We have shown that the multidimensional quadrature problem is 
solved exactly for any polynomial of degree two or less when using quadratic resampling. 
When tested on typical multidimensional asset pricing problems, quadratic resampling yields 
significant error reduction. The error reduction factors ranged from 1.2 to 1000 in our 
experiments. The typical error reduction for problems in up to 10 dimensions was about 3, 
corresponding to a speedup factor of roughly one order of magnitude. The method is readily 
applicable to any European asset pricing problem involving an arbitrary number of variables. 
We have successfully implemented the method for problems with up to 100 dimensions. 
Besides its applications to security pricing, the method could be applied to any problem where 
the computation of high-dimensional integrals is required. Examples of such problems can 
be found for example in failure rate analysis of mechanical and electronic systems, and in 
the simulation of complex physical phenomena such as those occurring in nuclear physics, 
mechanics, and fluid dynamics. 

We feel that the method presented in this paper and the experimental results obtained with 
it make it possible to realistically envision the use of multidimensional stochastic models for 
practical real-world quantitative risk management problems. This capability of computing 
the joint influences of several tens of risk factors such as interest rates of various terms in 
different currencies, equity of various kinds, and any other relevant economic variables, may 
dramatically increase the competitive advantage of quantitative methods over more traditional 
analysis techniques. An area of particular interest is that of estimating precisely the economic 
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value (net present value) of firms such as banks and insurance companies from a quantitative 
analysis of their balance sheet. However, this numerical pricing method does not generalize 
simply to American asset pricing problems, which cannot be reduced to multidimensional 
quadrature problems. Extensions of this approach for American asset pricing are left for future 
research. 
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